/* Routine to pick a massive star from the Kroupa (2001) IMF, which is 
   dN / dlog m \propto 
    m^0.7 (0.01 < m < 0.08)
    m^-0.8 (0.08 < m < 0.5)
    m^-1.7 (0.5 < m < 1)
    m^-1.3 (1 < m < 120)
*/ 
 
/* Probability of being in each interval of the Kroupa IMF */
  //#define wgt1 0.4967
  //#define wgt2 0.4360
  //#define wgt3 0.0426
  //#define wgt4 0.0247



extern double cluster_mass;

double igimf_imf(gsl_rng *rand) {
  
  double wgt1=0.;       //We exclude the low mass since there are no models for those
  double wgt2=0.866;	//we need to dynamically alter these weights
  double wgt3=0.0846;	//based on the maximum stellar mass...
  double wgt4=0.049;
  
  double x, mmin, mmax, gamma;
  double max_stellar_mass_igimf, delta_wgt;
  double pow00, x1, x2, x3, x4, x5,x6;
 // if (cluster_mass < 1E3)
//	max_stellar_mass_igimf=pow(10., 0.5792 * log10(cluster_mass) - 0.1762);
//  else 	max_stellar_mass_igimf=pow(10., 0.4867 * log10(cluster_mass) - 0.0209);
//taken from Carsten Weidner Email
  x1=log10(cluster_mass);
  x2=x1*x1;
  x3=x2*x1;
  x4=x3*x1;
  x5=x4*x1;
  x6=x5*x1;
  pow00= 1.448914324632-2.521767207421*x1+  2.054906559817*x2 -0.615883506965*x3+  0.089748388302*x4 -0.006426669425*x5+  0.000181597847*x6;
  max_stellar_mass_igimf=pow(10., pow00);

  if(max_stellar_mass_igimf > 120){
    max_stellar_mass_igimf=120;
  }
  else
    {
      //calculate how much less likely it is to fall in
      //the highest mass bin due to the fact that
      //we have a maximum mass that isn't 120  
      delta_wgt= pow(max_stellar_mass_igimf, (double) -1.3)/1.3;
      delta_wgt-= pow((double) 120, (double) -1.3)/1.3;
      
      delta_wgt/= pow((double) 1, (double) -1.3)/1.3 - 
	pow((double) 120, (double) -1.3)/1.3;
      
      delta_wgt*=wgt4;
      
      wgt2*=1./(1.-delta_wgt);
      wgt3*=1./(1.-delta_wgt);
      wgt4=wgt4-delta_wgt;
      wgt4*=1./(1.-delta_wgt);
    }
  
  //First decide which power-law interval to draw from. Note we use
  //gamma 1 less than the exponent in Kroupa (2001), because we want
  //to draw from dN / dm = (dN / dlog m) / m
 
  
  x = gsl_rng_uniform (rand);
  
  if (x <= wgt1) {
    mmin = 0.01;
    mmax = 0.08;
    gamma = -0.3;
  } else if (x <= wgt1+wgt2) {
    mmin = 0.08;
    mmax = 0.5;
    gamma = -1.8;
  } else if (x <= wgt1+wgt2+wgt3) {
    mmin = 0.5;
    mmax = 1.0;
    gamma = -2.7;
  } else {
    mmin = 1.0;
    mmax = max_stellar_mass_igimf;
    gamma = -2.3;
  }
  
  
  
  // Now draw from a power-law in m and return 
  return(ranpowerlaw(mmin, mmax, gamma, rand));
}
